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Abstract The effects implied for the structure of com- 
pact objects by the modification of General Relativity 
produced by the generalization of the Lagrangian den- 
sity to the form f(R) = R + aR 2 , where R is the Ricci 
curvature scalar, have been recently explored. It seems 
likely that this sguared-gravity may allow heavier Neu- 
tron Stars (NSs) than GR. In addition, these objects 
can be useful to constrain free parameters of modificd- 
gravity theories. The differences between alternative 
gravity theories is enhanced in the strong gravitational 
regime. In this regime, because of the complexity of the 
field equations, perturbative methods become a good 
choice to treat the problem. Following previous works 
in the field, we performed a numerical integration of the 
structure equations that describe NSs in /(i?)-gravity, 
recovering their mass-radius relations, but focusing on 
particular features that arise from this approach in the 
profiles of the NS interior. 

We show that these profiles run in correlation with 
the second-order derivative of the analytic approxima- 
tion to the Equation of State (EoS), which leads to re- 
gions where the enclosed mass decreases with the radius 
in a counter-intuitive way. We reproduce all computa- 
tions with a simple polytropic EoS to separate zcroth- 
order modified gravity effects. 
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1 Introduction 



neutron stars • gener- 



Current cosmological observations interpreted in the 
standard cosmological model require the presence of a 
non-standard matter content in order to explain the ac- 
celerated expansion of the Universe [H, [l7 , 18 , 21, 23| . 
Along the last decade alternative cosmological models 
have been developed to reinterpret these data with- 
out involving any unknown, doubtful component of the 
energy- matter tensor (see, for instance, |4|). The ap- 
pearance of Extended Theories of Gravity (ETGs) was 
strongly stimulated by the possibilities they might pro- 
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vide in this context [5 

ETGs are based on corrections and generalizations 
of Einstein's General Relativity (GR) theory. We fo- 
cus on a particular class, called /(-R)-gravity theories, 
which are based on a modification of the Einstein- Hilbcrt 
action: the usual Lagrangian density is generalized re- 
placing the Ricci curvature scalar R by a function of 
it or high-order invariants of the curvature tensor ,such 



as R 2 , R^W, R^ a f}R^ uaf> , RDR, RD k R (see [5( for 
a complete review). 

Several of these models are successfully constructed 
to satisfy the current Solar System and laboratory tests 
[H,[l6, 1!^ 28| . In particular, the simplest choice f{R) = 
R + aR 2 , also called "i?-squared" gravity, has been fur- 
ther studied as the basis for a viable alternative cos- 
mological model, that can lead to the accelerated ex- 
pansion of the Universe and is well consistent with the 
temperature anisotropics observed in Cosmological Mi- 
crowave Background But in contrast to gravity in 
the weak-field regime, which has been subject to nu- 



2 



M. Orellana et al. 



merous experimental tests, gravity in the strong-field 
regime is largely unconstrained by observations [e.g.Q. 

However, other authors, making a more detailed model 
of the structure of compact stars, obtained a set of mod- 
ified Tolman-Oppcnhcimcr-VolkofF (TOV) equations that 
describe a spherically symmetric mass distribution, un- 
der hydrostatic equilibrium, in simple /(i?)-gravity, and 
derived solutions that reproduce the correct behaviour 
at the weak gravity limit [e.g.0, in the context of scalar- 
tensor theories of gravity]. Furthermore, working with 
a perturbative approach to solve the field equations, 
Cooney et al. [3] and Santos [25| have found mass-radius 
relations for compact stars using a polytropic Equation 
of State (EoS). More recently Arapoglu et al. [2| and 



Deliduman et al. 10( applied the same approach using 



a set of realistic EoS for Neutron Stars (NSs) in the 
i?-squared and R^R^ gravities, respectively. 

The mass-radius relations obtained by Arapoglu et al. 
0] and Deliduman et al. [ltj indicate that such f(R) 
models can accommodate NSs up to masses larger than 
the currently observed ones, which are at most M max = 
1.97±O.O4M for PSR J1614-2230 The i?-squared 
gravity introduces a new parameter in the model through 
the value of a, the coefficient of R 2 . The freedom in 
the choice of this parameter allows some EoS, which 
are excluded within the framework of GR, to be rec- 
onciled with the observations. Motivated by these re- 
sults, we investigate in detail the structure of NSs un- 
der this model for two different EoS. One of them is a 
polytropic approximation that we use here to separate 
zeroth-order modified-gravitational effects, whereas the 
other provides a realistic representation of nuclear mat- 
ter at high densities. 

The paper is organized as follows. In Section 2, we 
obtain the modified TOV equations following the per- 
turbative approach to solve the field equations. In Sec- 
tion 3, we present the EoS used to integrate the equa- 
tions and we briefly describe the numerical methods of 
resolution. In Section 4, we present our results for the 
mass-radius relations, focusing on the behaviour of the 
profiles obtained for the NS interior. Final remarks are 
shown in Section 5. 



2 Structure equations in Ji-squared gravity 

The equations that define the structure of a NS in GR 
are deduced proposing the static, spherically symmetric 
line element, ds, to be 



equations, so-called TOV because of the pioneer work 
of and [2^], can be expressed as: 

dm(r) 



47rr p(r) 
dr 

dp(r) c 2 p(r)+p(r) 



ds 2 



Then, considering the Einstein equations for an ideal 
energy-momentum tensor T£ = diag{— pc 2 ,p,p,p}, these 



dr 2Gm(r) — c 2 
-2A(r) _ 1 2Gm{r) 



1 



d<P{r) 



G ( — r p{r) H 



dp 



(2) 
(3) 

(4) 
(5) 

dr (c 2 p(r) + p(r)) dr ' 

where m(r) is the total relativistic mass enclosed in 
a sphere of radius r. The functions p(r) and p(r) are, 
respectively, the mass-energy density and the pressure 
at this radius. We explicitly keep the gravitational con- 
stant, G, and the velocity of light, c, since the quantities 
are considered with their full-dimensions for integration 
purposes. Giving an explicit relation between p and p, 
the so-called EoS, the TOV equations can be solved as- 
suming a value for the central density, p(r = 0) = p Cl 
and integrating the system to the pressure vanishes, 
p(r = = 0. Here i?* is the radius of the star and 
= m(Rn) the stellar mass. Every p c generates a cou- 
ple of values M* and R+. Then, varying this parameter, 
a mass-radius (M* — i?*) relation is defined for every 
EoS. 

The modified TOV equations can be obtained from 
the gravitational field equations. Adding the new term 
to the Hilbert-Einstcin action, we have: 



S 



16nG 



d A x, 



^(R + aR 2 ) + Scatter, 



(0) 



where g is the determinant of the metric g^ v . 

Working in the metric formalism, the variation of 
the action with respect to the metric yields fourth-order 
differential equations of g^. This poses an enormous 
obstacle to solve the problem thoroughly. For this rea- 
son, we adopt the perturbative approachpresented by 
Cooney et al. Q and Arapoglu et al. Rewriting 
f(R) = R + aR 2 = R{1 + /?), we consider the f(R) 
function as a perturbation of a GR background. Hence, 
the dimensionless quantity /3 = aR comprises the de- 
viation from GR and the perturbative method can be 
applied as long as |/3| ^ 1. Under this condition, we 
can work with equations of motion up to first order 
in (3 without imposing any constrain at the level of the 
action and ensuring the nature of the variational princi- 
ple Neglecting terms with 0(/3 2 ) or higher, the field 
equations reduce, for this particular choice of f(R), to 



-e 2 *^c 2 dt 2 +e 2A ^dr 2 +r 2 (d6 2 +sm 2 8d V 2 ). (1) G u + a 



2R [ R^ v — -Rg^v 



2 [g^UR - V M V,i?) 



8ttG 
T 

4 ^1,, 



(7) 
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where G M „ = — \Rg^v is the Einstein tensor. In 
the limit a — > 0, /3 — > 0, and the field equations of 
GR are recovered [see e.g. [5[. 

We also assume a static and spherically symmetric 
line element, given by equation ((T|). The perturbativc 
approach allows to expand the functions present in the 
metric into a leading term (unperturbed) , denoted with 
subscript 0, plus a corrective one, denoted with sub- 
script 1, that is of first order in j3: A = Aq + f3A\ and 
<P = <Po+(3<I>i. The hydrodynamic quantities are also de- 
fined perturbatively: p = po + (3pi and p = po + j3pi 0. 
Hence, new restrictions are imposed over the value of 
/3 by the constraints (3<Pi < <P , $A\ < A Q , j3pi < po, 
and (3pi <C po- From now on, we use the prime for radial 
derivatives. 

Following Q and Q, we define the mass assuming 
that the solution for the metric has the same form as 
the exterior Schwarzschlid solution in GR, i.e. 



_ 2A{r) = 1 2GM. for r > R*. 



For the interior solution, 



-2A{r) 



1 



2Gm{r) 



for r < R+, 



(S) 



(9) 



Mass-derivative terms are indicated, giving: dm/dr = 
inr 2 p — f3[A+ B + C] . We study each term contribution 
below. 

Exact equations (|T0|) and (fTTT) arc impractical be- 
cause /3 involves other derivatives through _R, which 
makes (3 r-dependent in a complicated wajQ. Thus, we 
integrate equations (JTUJ) and (fTTj) assuming that /3 is 
well approximated by /3 ~ /3 = aRo, where Rq is the 
Ricci scalar locally calculated in GR and a is a constant 
parameter with squared distance units, compatible with 
other authors approach [for instance [2|. Note that 



Ro 



8Gtt 



{p c 2 - 3p ), 



(12) 



and then, contrary to the GR case, the derivatives of 
the EoS, dp /dp and d 2 p/dp 2 , also enter into (fit)]) and 
(fTTj) . through R' and R„. 

In /(i?)-gravity the weight of the perturbation is ad- 
justed by the value of the a parameter. In our work, we 
restrict ourselves to the constraints reported by Santos 



24j, and references therein, which points to 10 cm < 



a/3 < 10 cm , based on astronomical observations 
and nuclear experiments in terrestrial laboratories. 



where m also admits a perturbative expansion m = 
mo + /3mi, with mo the zeroth-order mass that is ob- 
tained integrating @. 

With this considerations, and taking into account 
that po and po satisfy Einstein's equations, the derived 
modified TOV equations are: 



A/2 



dm 2 
— — = Anr p - 
dr 



2/3 



4?rr po - —r R 



2irp r J - —r + ^-m I ^ - r I — r -m ]^ 



G 



Ro 



2G 



Ro 



B/2 



-C/2 



(10) 



dp 
dr 



c p + p 
2Gm — c 2 r 

„2 



G 



4tt 



-r p + 



-2/3 



4tt 



r 2 po+ 



8G 



r 2 R 



2- 



-Por 



G 



-mo 



R 'o 
R» 



(11) 



Note that 2/3 ... indicate the first order correction in 

/3 into the gradients dm/dr and dp/dr, respectively. In 
order to work up to first order in /3, terms between 
the brackets have been evaluated at order zero. It is 
important to note that dp/dr does not depend on i? '. 



3 Equations of state and numerical methods 

To solve the system of equations given by (flU)) and pT|) 
is necessary to bring a relation between the pressure 
p and the density p or the energy density e, the so- 
called EoS. The EoS contains the information of the 
behaviour of matter inside NSs through several orders 
of magnitude in density. Because the properties of the 
matter at the highest densities in the central region of 
NSs are not well understood, different EoS are proposed 
and then constrained with observations of masses and 
radii of actual NSs. 

An analytical representation of the EoS is required 
for solving the structure of NSs in modified theories of 
gravity, where hydrostatic equilibrium equations are of 
fourth-order. In such cases the usual interpolation tech- 
nique fails to accurately represent high-order deriva- 
tives 13| . Analytical representations of several EoS have 



been obtained through a consistent procedure by Haensel 
and Potekhin 14 1, who calculated the best- fit coeffi- 
cients of a polynomial expansion both in the crust and 
core density regimes. However, it must be emphasized 
that these analytical EoS are approximations obtained 

1 The Ricci scalar in terms of the functions of the metric 
CE) is: 



R = 



2 e- 2A [r 2 {<$'A' + (<2>') 2 + <2>"} + 2r (<?>' -A 1 )- e 2A + l] 
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Fig. 1 Mass-radius (M* — i?*) relations for the two selected EoS: SLY and POLY (left and right, respectively), considering 
seven values for the a parameter, which are indicated above in km 2 units. All the curves correspond to values of central density, 
p c , in the range 10 14 6 — 10 15 9 gr cm -3 . The crosses indicate the maximum mass for each curve, assuming a necessary condition 
for equilibrium: dM/dp c > 0. 



by fitting only the zeroth-order relation between p and 
p, because it is the relation needed to calculate the 
structure of NSs in GR. Thus, special care should be 
taken if high-order derivatives of these expressions are 
used during the calculation, as in the case we are inter- 
ested in here (i.e. dp/ dp and d 2 p/dp 2 ). 

Taking this into account, and in order to compare 
our results with those already published in the litera- 
ture, we calculate mass-radius (Af* — i?*) relations con- 
sidering two different EoS: SLY [HEl] and POLY [26 1. 

The first one is a realistic EoS that properly repre- 
sents the behaviour of nuclear matter at high density. 
Its analytic parametrization is given by 



ai +a 2 £ + a 3 £ 3 

C = r- ~ c Jo (as 4 - a 6 )) 

1 + 04 £ 

+(a 7 + a 8 /o (09(010 - 0) 

+ (on + O12O /o(oi3(oi4 - 0) 

+(a i5 + aieO /o(ai7(ai8 - 0) > 



where 



C = log(p/gcm 3 ), C = log(P/dyncm" 



fo(x) 



1' 



(13) 



(14) 



(15) 



and the coefficients a, are tabulated [lJJ. This is the 
same expression used by and we use it here to test 



our results. The second EoS is a simpler polytropic ap- 
proximation given by 



C = 2^ + 5.29355 



(16) 



Despite the latter is not a realistic EoS that thoroughly 
represents NSs, it is a toy model that allows to study 
zeroth-order modified gravity effects, separating them 
from more tricky effects arising in the case of a real- 
istic EoS, with its complex analytical expression and 
from which the error propagating to the derivatives is 
out of our control. The precise value of the adiabatic 
index, r = dlogp/dlogp = d£/d£, is not relevant as r 
remains a derivable function. The reader is referred to 
22| for tighter constraints on r that point to a some- 
what larger value than the one adopted here. 

3.1 Numerical Method 

Solving the system of ordinary differential equations 
formed by the equations (fTU)) and (fTTj) implies their 
integration from the centre to the NS surface, using 
the chosen EoS. Once the solution is found, a couple 
of values for the mass, M*, and the radius, i?*, are es- 
tablished. In order to perform the integration, we use 
a numerical code based on a fourth-order Runge-Kutta 
method on the radial coordinate. For this coordinate 
we implement a variable step which is systematically 
shortened close to the NS surface, to account for rapid 
variations of the physical parameters in this region. 
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r [km] r [km] 

Fig. 2 Profiles of the internal structure of NSs for two extreme cases of low and high central densities, p c = 10 14,6 and 
10 15,4 gr cm -3 , and for three different values of the a parameter (+0.2, 0.0 and -0.2 km 2 ), where a. = 0.0 corresponds to GR 
case. On the left (right) panel profiles corresponding to the SLY (POLY) EoS are shown. A zoom-in of the mass profile close 
to the NS surface is shown in Figure [3] At low central density values the effect on the integrated mass can still represent a 
deviation as much as < 10% from the GR mass. 




Fig. 3 Zoom-in of the the mass and density profiles close of the surface of the NS for the SLY EoS shown in Figure [2] for 
p c = 10 15 ' 4 and 10 14,6 gr cm -3 . For the value of the a parameter -0.2 km 2 (+0.2 km 2 ) the mass increases (decreases) roughly 
10% respect to de GR case (a = 0), for both low and high central densities. 



During the Rungc-Kutta loop, we also solve the dif- 
ferential equations corresponding to the metric compo- 
nents: g t t and g rr , which are involved in the criteria for 
the validity of the perturbative approach. In each step, 
we first integrate the TOV equations in the frame of 
GR to obtain zeroth-order values that then we use to 
calculate the first-order perturbative solution. 

We start the numerical integration from the centre 
with a given central density, p c , and we finish the inte- 
gration at the surface, defining the NS radius, R±, and 
mass, M*, when the density reaches p = 10 6 gr cm -3 . 
This density corresponds to the outer boundary of the 



NS crust and is the limit of validity for these kind of 
EoS, as they were conceived beginning with a model for 
nuclear matter at high densities. 

4 Results 

In Figure Q] we present the mass-radius (M* — R+) rela- 
tions obtained for SLY and POLY EoS (left and right 
panels, respectively), using seven values of the a pa- 
rameter between -0.2 and +0.2 km 2 and considering 
central densities, p c , ranging from 10 14 ' 6 gr cm' 3 to 
10 15 - 9 gr cm" 3 . Maximum masses achieved are indi- 
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cated by crosses in each curve, assuming a necessary 
condition for equilibrium: dM/dp c > 0. Our results for 
SLY EoS are in accordance with those previously pre- 
sented by [2| • Values of a < km 2 (a > km 2 ) can ac- 
commodate higher (lower) maximum masses than GR. 
In particular, POLY configurations are less sensitive to 
the value of a than those from SLY EoS. 

In Figure[5]we present the internal profiles found for 
the density, p{r), and mass, m(r), for our a extreme val- 
ues and the GR case. The density and the pressure fol- 
low rather usual (resembling GR) profiles, where both 
magnitudes monotonously decrease with radius. How- 
ever, particular effects of f{R) are reflected in the mass 
profiles for the realistic SLY EoS, and become more pro- 
nounced for the high-mass NSs (p c — 10 15 4 gr cm -3 ). 

These effects are evident close to the NS surface (at 
r ~ 10 km) where, in a narrow layer {Ar ~ 0.2 km), 
an unexpected (counter- intuitive) decrease in m(r) ap- 
pears before (a > km 2 ) or after (a < km 2 ) a dip 
(peak) in the profile. In Figure |3] we present a zoom of 
the mass and density profiles (left and right panels, re- 
spectively) close to the NS surface, obtained using the 
SLY EoS for logp c [gr cm~ 3 ] = 15.4 and 14.6. Besides 
the modification of M* is higher for high central density 
stars, the relative change in the total mass with respect 
to the GR case is roughly 10% in both high/low central 
density cases. 

In the frame of GR, a decreasing mass profile could 
only be accomplished by a fluid of negative density, 
because dm/dr = 47rr 2 p(r). However, in /(i?)-gravity 
these profiles can be explained as a consequence of the 
modified geometry. In contrast, no such features are 
present in the profiles obtained with the POLY EoS for 
all the values considered for the parameters. 

In Figurc[5]wc present the profile of the ratio g^ r /g rr 
with the density through all the NS interior for the SLY 
EoS, considering the same densities of Figure We re- 
call that this ratio should be close to 1.0 as a necessary 
condition for the validity of the pcrtubative method. On 
the upper panel of this figure we also plot the first and 
second logarithmic derivatives of the SLY EoS. From 
the comparison of the trend, a strong coupling of the 
perturbative deviations of g rr with the second-order 
derivative of the EoS is evident. The modification in the 
metric radial component is mild, and only perceptible 
when the second-order derivative becomes important, 
strongly oscillating, in the 10 11 — 10 14 gr cm -3 density 
range. Such behaviour is not present for the POLY EoS, 
which logarithmic second derivative is null, maintaining 
9rr/drr — 1 all through the NS interior. It is important 
to note that the function <P(r) in the time component 
of the metric is actually reflecting the behaviour of the 
pressure, whereas A{r) is governed by the mass, and 
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Fig. 4 Top panel: Profile of the of the first (dotted line) and 
second (continuous line) logarithmic derivatives of the SLY 
EoS in the NS interior. Middle and Bottom panels: Profiles 
of the ratio between the radial component of the metric at or- 
der zero (GR), and at first order {g® r /g r r) for a = +0.2 and 
—0.2 km 2 , which should be close to 1.0 as a necessary condi- 
tion of the pertubative method. The perturbative deviations 
are closely related with the behaviour of the second-order 
derivative of the EoS. 



it is the mass but not the pressure the one requiring 
evaluations of Rq, and thus depending on high-order 
derivatives of the EoS. 

To further analyse the origin of the deviations from 
the GR we explore the contribution to dm/dr of the 
four terms involved in equation (|10j) . which we call: 
A-Kr 2 p (GR) and A, B, C (perturbative terms). In Fig- 
ure[5]we present each term contribution when a = +0.2 
km 2 for SLY (upper-left panel) and POLY (upper-right 
panel) EoSs, in the case of a high central density star 
(logpc [gr cm -3 ] = 15.4). For the SLY EoS and for val- 
ues of radii r < 9.5 km, dm/dr is dominated by the GR 
term. Closer to the NS surface, for r > 9.5 km, the term 
C oc Rq/Rq becomes dominant, causing the fluctuation 
in the mass profile. On the contrary, for the POLY EoS, 
which derivatives are strictly bounded smooth func- 
tions, this counter-intuitive effect is not present and 
the trend of dm/dr is dominated by the GR term, with 
very small modifications due to the perturbative terms. 
In the lower panels of Figure we zoom-in the upper 



Structure of neutron stars in i?-squared gravity 



7 



SLY, a= +0.2, log p c =15.4 



POLY, a= +0.2, log p c =15.4 
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Fig. 5 Profile of the mass gradient, dm/dr (thick line), close to the surface of the NSs for SLY (left panel) and POLY (right 
panel) EoSs. Deviations of the mass profile from the GR case are much more importante for the realistic SLY EoS. Note that 
the dashed line indicates zeroth-order (GR) term and the continuous with plus signs line indicates the contribution of the 
C oc Rq /Rq term. Lower panels zoomed-in to show in detail the contribution of the minor perturbative terms, indicated in 
the legend. 



SLY, a= -0.2, log p c =15.4 




4xicr 



-4x10" 



SLY, a= +0.2, log p c =1 4,6 
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Fig. 6 Same as Figure [5] for other three SLY cases: log p c [gr cm -3 ] = 15.4 (left) and 14.6 (center and right). The central 
and right plots compare the effects of changing the sign of a, which is indicated in km 2 units. The same values, i.e. a = —0.2 
(left) and +0.2 (Fig[5j left pannel) are shown for a large central density. 



panels to show the behaviour of the minor perturbative 
terms. 

In Figure [HI we extend this analysis to compare the 
behaviour of the mass derivative in four SLY cases: 
logp c [gr cm~ 3 ] = 15.4 and 14.6 for a = +0.2 and 
-0.2 km 2 . For low central density stars, the effect of 
the second-order derivative (C term), is relatively less 
important than in the high central density case. The 
fluctuations in this term occur in the density range 
where the second order logarithmic derivatives of the 
EoS become relevant, which in this case corresponds to 



a wider radial band 7-10.5 km, as the density changes 
in smoother way than in the high central density case. 
In the latter the density drastically decays in a narrow 
and superficial range from 9.7 to 10.2 km. 



In all cases studied here, the other terms, namely A 
and B, which correspond to ~ Rq and 
tions, respectively, are orders of magnitude lower. 



R' contribu- 
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5 Conclusions 

With the aim to investigate whether /(i?)-theories are 
viable to describe astrophysical scenarios like NSs, we 
have studied the particular i?-squared case using both 
simplified and realistic EoS. We have followed the gen- 
eral steps presented by Cooney et al. Q and Arapoglu 
et al. .2] using a perturbative approach applied to solve 
the fourth order field equations. Concerning the mass- 
radius — R+) relations, we have obtained results 
consistent with former studies, finding that for the high- 
est absolute values admitted for the a parameter, f(R)- 
thcories can accommodate heavier NSs than GR for ev- 
ery EoS. In this sense, it is important to remark that 
there is no agreement on the maximum mass achiev- 
able by NSs before they collapse to black holes, based 
on the uncertainties present on the behaviour of nuclear 
matter at the highest densities through their EoS. This 
problem can not be split out from our lack of under- 
standing of gravity |30t | . 

Notwithstanding, our most notorious result concerns 
the details of the internal structure of NSs considering 
the largest acceptable value for the a parameter (i.e. 
the stronger perturbation allowed to GR by the actual 
constrains). We find that the behaviour of the metric, 
which in i?-squarcd gravity depends not only on the 
EoS, but also on its higher-order derivatives, leads to 
a region where the mass enclosed decreases with the 
radius. 

Despite the fact that, in the frame of GR, this ef- 
fect could only be explained by means of a negative- 
density fluid, in f(R) theories, it arises as a natural 
consequence of the coupled space-time geometry and 
matter content. Adopting a simple polytropic EoS, with 
strictly bounded higher-order derivatives, no anomalous 
behaviours of the internal profiles of the NS structure 
are observed, and the final effect of modified gravity is 
reduced. 

We emphasize that the uncertainties on the adopted 
EoS could have an enhanced effect on the solutions. 
Therefore we add a word of caution, as it remains un- 
clear whether the spikes in the mass profiles arise as a 
consequence of the analytical representation of the EoS, 
the perturbative approach or the geometry of squared- 
gravity. In any case, our results raise new questions 
in the topic. Further research is needed to disentangle 
these different possibilites. 

Finally, we suggest that a thorough study of the 
stability of the calculated structures would be very im- 
portant to ensure that these particular configurations 
can be realized under i?-squarcd gravity. Such a work, 
as well as the study of the astrophysical implications, 
is beyond the scopes of this paper. 
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